module obs_raob_mod

  use datetime
  use container
  use flogger
  use fiona
  use string
  use data_translators_utils_mod
  use latlon_process_mod
  use latlon_parallel_mod
  use raob_mod
  use raob_cimiss_xml_mod
  use const_mod
  use namelist_mod
  use latlon_process_mod
  use obs_operator_types_mod
  use obs_err_mod
  use efso_utils_mod

  implicit none

  private

  public obs_raob_type
  public raob_record_type

  type obs_raob_type
    integer :: obs_size       = 0
    integer :: max_man_level  = 0
    integer :: max_sigt_level = 0
    integer :: max_sigw_level = 0
    integer :: max_trop_level = 0
    type(hash_table_type), allocatable :: platforms
    type(linked_list_type), allocatable :: records                   ! Yo
    type(obs_grid_connector_type), allocatable :: connectors_man (:)
    type(obs_grid_connector_type), allocatable :: connectors_sigt(:)
    type(obs_grid_connector_type), allocatable :: connectors_sigw(:)
    type(obs_grid_connector_type), allocatable :: connectors_trop(:)
    type(raob_record_type), allocatable :: simulated_mean_records(:) ! H mean(X)
    type(raob_record_type), allocatable :: mean_simulated_records(:) ! mean(Yb)
    type(raob_record_type), allocatable :: simulated_records(:,:)    ! Yb or Yb'
    type(raob_record_type), allocatable :: innov_records(:)          ! Yo - H mean(X)
    type(obs_err_type), allocatable :: obs_err(:)
    integer, allocatable :: mask(:)
  contains
    procedure :: init => obs_raob_init
    procedure :: read => obs_raob_read
    procedure :: filter_by_latlon_and_time => obs_raob_filter_by_latlon_and_time
    procedure :: clear => obs_raob_clear
    procedure :: calc_mean_simulated_records => obs_raob_calc_mean_simulated_records
    procedure :: calc_perb_simulated_records => obs_raob_calc_perb_simulated_records
    procedure :: calc_innov => obs_raob_calc_innov
    procedure :: init_write => obs_raob_init_write
    procedure :: write => obs_raob_write
    final :: obs_raob_final
  end type obs_raob_type

  real, parameter :: raob_man_p(16) = [ &
    1000.0, 925.0, 850.0, 700.0, 500.0, 400.0, 300.0, 250.0, 200.0, 150.0, 100.0, 70.0, 50.0, 30.0, 20.0, 10.0 &
  ]

contains

  subroutine obs_raob_init(this)

    class(obs_raob_type), intent(inout) :: this

    call this%clear()

    allocate(this%platforms)
    allocate(this%records  )

    this%platforms = hash_table(chunk_size=100000)

  end subroutine obs_raob_init

  subroutine obs_raob_read(this, file_path, mute)

    class(obs_raob_type), intent(inout) :: this
    character(*), intent(in) :: file_path
    logical, intent(in) :: mute

    call raob_cimiss_xml_read(file_path, this%platforms, this%records, mute)

  end subroutine obs_raob_read

  subroutine obs_raob_filter_by_latlon_and_time(this, start_lon, end_lon, start_lat, end_lat, start_time, end_time)

    class(obs_raob_type), intent(inout) :: this
    real(r8), intent(in) :: start_lon
    real(r8), intent(in) :: end_lon
    real(r8), intent(in) :: start_lat
    real(r8), intent(in) :: end_lat
    type(datetime_type), intent(in) :: start_time
    type(datetime_type), intent(in) :: end_time

    type(linked_list_iterator_type) it
    integer n, i, j, k

    it = linked_list_iterator(this%records)
    do while (.not. it%ended())
      select type (record => it%value)
      type is (raob_record_type)
        if (record%station%lon < start_lon * deg .or. record%station%lon > end_lon * deg .or. &
            record%station%lat < start_lat * deg .or. record%station%lat > end_lat * deg .or. &
            record%time < start_time .or. record%time > end_time) then
          call this%records%remove_item(it%item)
        end if
      end select
      call it%next()
    end do

    if (parallel_ensemble .and. proc%is_root_proc()) then
      call log_notice('After filter there is ' // to_str(this%records%size) // ' raob records.')
    end if

    if (allocated(this%connectors_man)) deallocate(this%connectors_man)
    allocate(this%connectors_man(this%records%size))
    if (allocated(this%connectors_sigt)) deallocate(this%connectors_sigt)
    allocate(this%connectors_sigt(this%records%size))
    if (allocated(this%connectors_sigw)) deallocate(this%connectors_sigw)
    allocate(this%connectors_sigw(this%records%size))
    if (allocated(this%connectors_trop)) deallocate(this%connectors_trop)
    allocate(this%connectors_trop(this%records%size))
    if (allocated(this%simulated_mean_records)) deallocate(this%simulated_mean_records)
    allocate(this%simulated_mean_records(this%records%size))
    if (allocated(this%mean_simulated_records)) deallocate(this%mean_simulated_records)
    allocate(this%mean_simulated_records(this%records%size))
    if (allocated(this%simulated_records)) deallocate(this%simulated_records)
    allocate(this%simulated_records(this%records%size,ens_ibeg:ens_iend))
    if (allocated(this%innov_records)) deallocate(this%innov_records)
    allocate(this%innov_records(this%records%size))

    n = 0
    it = linked_list_iterator(this%records)
    do i = 1, this%records%size
      select type (record => it%value)
      type is (raob_record_type)
        n = n + record%man %num_level; this%max_man_level  = max(this%max_man_level , record%man %num_level)
        n = n + record%sigt%num_level; this%max_sigt_level = max(this%max_sigt_level, record%sigt%num_level)
        n = n + record%sigw%num_level; this%max_sigw_level = max(this%max_sigw_level, record%sigw%num_level)
        n = n + record%trop%num_level; this%max_trop_level = max(this%max_trop_level, record%trop%num_level)
      end select
      call it%next()
    end do
    this%obs_size = n

    ! Initialize observation error structure.
    allocate(this%obs_err(this%obs_size))
    n = 1
    it = linked_list_iterator(this%records)
    do i = 1, this%records%size
      select type (record => it%value)
      type is (raob_record_type)
        ! Mandatory levels
        do j = 1, record%man%num_level
          if (is_missing(record%man%p(j))) call log_error('RAOB mandatory pressure is missing!', __FILE__, __LINE__)
          do k = 1, size(raob_err_plev)
            if (raob_err_plev_bnds(1,k) <= record%man%p(j) .and. record%man%p(j) < raob_err_plev_bnds(2,k)) exit
          end do
          if (k > size(raob_err_plev)) call log_error('RAOB does not found observation error for level ' // to_str(int(record%man%p(j))) // 'hPa!')
          this%obs_err(n)%u = u_err_raob(k)
          this%obs_err(n)%v = v_err_raob(k)
          this%obs_err(n)%t = t_err_raob(k)
          this%obs_err(n)%q = q_err_raob(k)
          this%obs_err(n)%p = p_err_raob(k)
          n = n + 1
        end do
        ! Significant temperature levels
        do j = 1, record%sigt%num_level
          if (is_missing(record%sigt%p(j))) call log_error('RAOB significant temperature pressure is missing!', __FILE__, __LINE__)
          do k = 1, size(raob_err_plev)
            if (raob_err_plev_bnds(1,k) <= record%sigt%p(j) .and. record%sigt%p(j) < raob_err_plev_bnds(2,k)) exit
          end do
          if (k > size(raob_err_plev)) call log_error('RAOB does not found observation error for level ' // to_str(int(record%sigt%p(j))) // 'hPa!')
          this%obs_err(n)%u = u_err_raob(k)
          this%obs_err(n)%v = v_err_raob(k)
          this%obs_err(n)%t = t_err_raob(k)
          this%obs_err(n)%q = q_err_raob(k)
          this%obs_err(n)%p = p_err_raob(k)
          n = n + 1
        end do
        ! Significant wind levels
        do j = 1, record%sigw%num_level
          if (is_missing(record%sigw%p(j))) call log_error('RAOB significant wind pressure is missing!', __FILE__, __LINE__)
          do k = 1, size(raob_err_plev)
            if (raob_err_plev_bnds(1,k) <= record%sigw%p(j) .and. record%sigw%p(j) < raob_err_plev_bnds(2,k)) exit
          end do
          if (k > size(raob_err_plev)) call log_error('RAOB does not found observation error for level ' // to_str(int(record%sigw%p(j))) // 'hPa!')
          this%obs_err(n)%u = u_err_raob(k)
          this%obs_err(n)%v = v_err_raob(k)
          this%obs_err(n)%t = t_err_raob(k)
          this%obs_err(n)%q = q_err_raob(k)
          this%obs_err(n)%p = p_err_raob(k)
          n = n + 1
        end do
        ! Tropopause levels
        do j = 1, record%trop%num_level
          if (is_missing(record%trop%p(j))) call log_error('RAOB tropopause pressure is missing!', __FILE__, __LINE__)
          do k = 1, size(raob_err_plev)
            if (raob_err_plev_bnds(1,k) <= record%trop%p(j) .and. record%trop%p(j) < raob_err_plev_bnds(2,k)) exit
          end do
          if (k > size(raob_err_plev)) call log_error('RAOB does not found observation error for level ' // to_str(int(record%trop%p(j))) // 'hPa!')
          this%obs_err(n)%u = u_err_raob(k)
          this%obs_err(n)%v = v_err_raob(k)
          this%obs_err(n)%t = t_err_raob(k)
          this%obs_err(n)%q = q_err_raob(k)
          this%obs_err(n)%p = p_err_raob(k)
          n = n + 1
        end do
      end select
      call it%next()
    end do

    ! Set mask for missing values.
    allocate(this%mask(this%obs_size * 4))
    it = linked_list_iterator(this%records)
    i = 1
    do while (.not. it%ended())
      select type (record => it%value)
      type is (raob_record_type)
        ! Mandatory levels
        do k = 1, record%man%num_level
          this%mask(i) = merge(0, 1, is_missing(record%man %ua(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%man %va(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%man %ta(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%man %sh(k))); i = i + 1
        end do
        ! Significant temperature levels
        do k = 1, record%sigt%num_level
          this%mask(i) = merge(0, 1, is_missing(record%sigt%ua(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%sigt%va(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%sigt%ta(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%sigt%sh(k))); i = i + 1
        end do
        ! Significant wind levels
        do k = 1, record%sigw%num_level
          this%mask(i) = merge(0, 1, is_missing(record%sigw%ua(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%sigw%va(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%sigw%ta(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%sigw%sh(k))); i = i + 1
        end do
        ! Tropopause levels
        do k = 1, record%trop%num_level
          this%mask(i) = merge(0, 1, is_missing(record%trop%ua(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%trop%va(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%trop%ta(k))); i = i + 1
          this%mask(i) = merge(0, 1, is_missing(record%trop%sh(k))); i = i + 1
        end do
      end select
      call it%next()
    end do

  end subroutine obs_raob_filter_by_latlon_and_time

  subroutine obs_raob_calc_mean_simulated_records(this)

    class(obs_raob_type), intent(inout) :: this

    real(r8), allocatable :: tmp(:,:)
    integer i, iens

    do i = 1, this%records%size
      this%mean_simulated_records(i)%man %ua = 0
      this%mean_simulated_records(i)%man %va = 0
      this%mean_simulated_records(i)%man %ta = 0
      this%mean_simulated_records(i)%man %sh = 0
      this%mean_simulated_records(i)%sigt%ua = 0
      this%mean_simulated_records(i)%sigt%va = 0
      this%mean_simulated_records(i)%sigt%ta = 0
      this%mean_simulated_records(i)%sigt%sh = 0
      this%mean_simulated_records(i)%sigw%ua = 0
      this%mean_simulated_records(i)%sigw%va = 0
      this%mean_simulated_records(i)%sigw%ta = 0
      this%mean_simulated_records(i)%sigw%sh = 0
      this%mean_simulated_records(i)%trop%ua = 0
      this%mean_simulated_records(i)%trop%va = 0
      this%mean_simulated_records(i)%trop%ta = 0
      this%mean_simulated_records(i)%trop%sh = 0
    end do
    do iens = ens_ibeg, ens_iend
      do i = 1, this%records%size
        ! Mandatory levels
        this%mean_simulated_records(i)%man %p  = this%mean_simulated_records(i)%man%p
        this%mean_simulated_records(i)%man %ua = add(this%simulated_records(i,iens)%man %ua, this%mean_simulated_records(i)%man %ua)
        this%mean_simulated_records(i)%man %va = add(this%simulated_records(i,iens)%man %va, this%mean_simulated_records(i)%man %va)
        this%mean_simulated_records(i)%man %ta = add(this%simulated_records(i,iens)%man %ta, this%mean_simulated_records(i)%man %ta)
        this%mean_simulated_records(i)%man %sh = add(this%simulated_records(i,iens)%man %sh, this%mean_simulated_records(i)%man %sh)
        ! Significant temperature levels
        this%mean_simulated_records(i)%sigt%p  = this%mean_simulated_records(i)%sigt%p
        this%mean_simulated_records(i)%sigt%ua = add(this%simulated_records(i,iens)%sigt%ua, this%mean_simulated_records(i)%sigt%ua)
        this%mean_simulated_records(i)%sigt%va = add(this%simulated_records(i,iens)%sigt%va, this%mean_simulated_records(i)%sigt%va)
        this%mean_simulated_records(i)%sigt%ta = add(this%simulated_records(i,iens)%sigt%ta, this%mean_simulated_records(i)%sigt%ta)
        this%mean_simulated_records(i)%sigt%sh = add(this%simulated_records(i,iens)%sigt%sh, this%mean_simulated_records(i)%sigt%sh)
        ! Significant wind levels
        this%mean_simulated_records(i)%sigw%p  = this%mean_simulated_records(i)%sigw%p
        this%mean_simulated_records(i)%sigw%ua = add(this%simulated_records(i,iens)%sigw%ua, this%mean_simulated_records(i)%sigw%ua)
        this%mean_simulated_records(i)%sigw%va = add(this%simulated_records(i,iens)%sigw%va, this%mean_simulated_records(i)%sigw%va)
        this%mean_simulated_records(i)%sigw%ta = add(this%simulated_records(i,iens)%sigw%ta, this%mean_simulated_records(i)%sigw%ta)
        this%mean_simulated_records(i)%sigw%sh = add(this%simulated_records(i,iens)%sigw%sh, this%mean_simulated_records(i)%sigw%sh)
        ! Tropopause levels
        this%mean_simulated_records(i)%trop%p  = this%mean_simulated_records(i)%trop%p
        this%mean_simulated_records(i)%trop%ua = add(this%simulated_records(i,iens)%trop%ua, this%mean_simulated_records(i)%trop%ua)
        this%mean_simulated_records(i)%trop%va = add(this%simulated_records(i,iens)%trop%va, this%mean_simulated_records(i)%trop%va)
        this%mean_simulated_records(i)%trop%ta = add(this%simulated_records(i,iens)%trop%ta, this%mean_simulated_records(i)%trop%ta)
        this%mean_simulated_records(i)%trop%sh = add(this%simulated_records(i,iens)%trop%sh, this%mean_simulated_records(i)%trop%sh)
      end do
    end do

    ! Communicate with other processes.
    ! Mandatory levels
    allocate(tmp(this%max_man_level,this%records%size))
    ! - u
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%man%num_level,i) = this%mean_simulated_records(i)%man%ua
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%man%ua = tmp(:this%mean_simulated_records(i)%man%num_level,i) / ensemble_size
    end do
    ! - v
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%man%num_level,i) = this%mean_simulated_records(i)%man%va
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%man%va = tmp(:this%mean_simulated_records(i)%man%num_level,i) / ensemble_size
    end do
    ! - t
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%man%num_level,i) = this%mean_simulated_records(i)%man%ta
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%man%ta = tmp(:this%mean_simulated_records(i)%man%num_level,i) / ensemble_size
    end do
    ! - q
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%man%num_level,i) = this%mean_simulated_records(i)%man%sh
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%man%sh = tmp(:this%mean_simulated_records(i)%man%num_level,i) / ensemble_size
    end do
    deallocate(tmp)
    ! Significant temperature levels
    allocate(tmp(this%max_sigt_level,this%records%size))
    ! - u
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%sigt%num_level,i) = this%mean_simulated_records(i)%sigt%ua
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%sigt%ua = tmp(:this%mean_simulated_records(i)%sigt%num_level,i) / ensemble_size
    end do
    ! - v
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%sigt%num_level,i) = this%mean_simulated_records(i)%sigt%va
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%sigt%va = tmp(:this%mean_simulated_records(i)%sigt%num_level,i) / ensemble_size
    end do
    ! - t
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%sigt%num_level,i) = this%mean_simulated_records(i)%sigt%ta
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%sigt%ta = tmp(:this%mean_simulated_records(i)%sigt%num_level,i) / ensemble_size
    end do
    ! - q
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%sigt%num_level,i) = this%mean_simulated_records(i)%sigt%sh
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%sigt%sh = tmp(:this%mean_simulated_records(i)%sigt%num_level,i) / ensemble_size
    end do
    deallocate(tmp)
    ! Significant wind levels
    allocate(tmp(this%max_sigw_level,this%records%size))
    ! - u
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%sigw%num_level,i) = this%mean_simulated_records(i)%sigw%ua
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%sigw%ua = tmp(:this%mean_simulated_records(i)%sigw%num_level,i) / ensemble_size
    end do
    ! - v
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%sigw%num_level,i) = this%mean_simulated_records(i)%sigw%va
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%sigw%va = tmp(:this%mean_simulated_records(i)%sigw%num_level,i) / ensemble_size
    end do
    ! - t
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%sigw%num_level,i) = this%mean_simulated_records(i)%sigw%ta
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%sigw%ta = tmp(:this%mean_simulated_records(i)%sigw%num_level,i) / ensemble_size
    end do
    ! - q
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%sigw%num_level,i) = this%mean_simulated_records(i)%sigw%sh
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%sigw%sh = tmp(:this%mean_simulated_records(i)%sigw%num_level,i) / ensemble_size
    end do
    deallocate(tmp)
    ! Tropopause levels
    allocate(tmp(this%max_trop_level,this%records%size))
    ! - u
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%trop%num_level,i) = this%mean_simulated_records(i)%trop%ua
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%trop%ua = tmp(:this%mean_simulated_records(i)%trop%num_level,i) / ensemble_size
    end do
    ! - v
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%trop%num_level,i) = this%mean_simulated_records(i)%trop%va
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%trop%va = tmp(:this%mean_simulated_records(i)%trop%num_level,i) / ensemble_size
    end do
    ! - t
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%trop%num_level,i) = this%mean_simulated_records(i)%trop%ta
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%trop%ta = tmp(:this%mean_simulated_records(i)%trop%num_level,i) / ensemble_size
    end do
    ! - q
    do i = 1, this%records%size
      tmp(:this%mean_simulated_records(i)%trop%num_level,i) = this%mean_simulated_records(i)%trop%sh
    end do
    call global_sum(proc%comm, tmp)
    do i = 1, this%records%size
      this%mean_simulated_records(i)%trop%sh = tmp(:this%mean_simulated_records(i)%trop%num_level,i) / ensemble_size
    end do
    deallocate(tmp)

  end subroutine obs_raob_calc_mean_simulated_records

  subroutine obs_raob_calc_perb_simulated_records(this)

    class(obs_raob_type), intent(inout) :: this

    integer i, iens

    ! Perturbation = Y - mean(Y)
    ! Mandatory levels
    do iens = ens_ibeg, ens_iend
      do i = 1, this%records%size
        this%simulated_records(i,iens)%man%p  = this%mean_simulated_records(i)%man%p
        this%simulated_records(i,iens)%man%ua = this%simulated_records(i,iens)%man%ua - this%mean_simulated_records(i)%man%ua
        this%simulated_records(i,iens)%man%va = this%simulated_records(i,iens)%man%va - this%mean_simulated_records(i)%man%va
        this%simulated_records(i,iens)%man%ta = this%simulated_records(i,iens)%man%ta - this%mean_simulated_records(i)%man%ta
        this%simulated_records(i,iens)%man%sh = this%simulated_records(i,iens)%man%sh - this%mean_simulated_records(i)%man%sh
      end do
    end do
    ! Significant temperature levels
    do iens = ens_ibeg, ens_iend
      do i = 1, this%records%size
        this%simulated_records(i,iens)%sigt%p  = this%mean_simulated_records(i)%sigt%p
        this%simulated_records(i,iens)%sigt%ua = this%simulated_records(i,iens)%sigt%ua - this%mean_simulated_records(i)%sigt%ua
        this%simulated_records(i,iens)%sigt%va = this%simulated_records(i,iens)%sigt%va - this%mean_simulated_records(i)%sigt%va
        this%simulated_records(i,iens)%sigt%ta = this%simulated_records(i,iens)%sigt%ta - this%mean_simulated_records(i)%sigt%ta
        this%simulated_records(i,iens)%sigt%sh = this%simulated_records(i,iens)%sigt%sh - this%mean_simulated_records(i)%sigt%sh
      end do
    end do
    ! Significant wind levels
    do iens = ens_ibeg, ens_iend
      do i = 1, this%records%size
        this%simulated_records(i,iens)%sigw%p  = this%mean_simulated_records(i)%sigw%p
        this%simulated_records(i,iens)%sigw%ua = this%simulated_records(i,iens)%sigw%ua - this%mean_simulated_records(i)%sigw%ua
        this%simulated_records(i,iens)%sigw%va = this%simulated_records(i,iens)%sigw%va - this%mean_simulated_records(i)%sigw%va
        this%simulated_records(i,iens)%sigw%ta = this%simulated_records(i,iens)%sigw%ta - this%mean_simulated_records(i)%sigw%ta
        this%simulated_records(i,iens)%sigw%sh = this%simulated_records(i,iens)%sigw%sh - this%mean_simulated_records(i)%sigw%sh
      end do
    end do
    ! Tropopause levels
    do iens = ens_ibeg, ens_iend
      do i = 1, this%records%size
        this%simulated_records(i,iens)%trop%p  = this%mean_simulated_records(i)%trop%p
        this%simulated_records(i,iens)%trop%ua = this%simulated_records(i,iens)%trop%ua - this%mean_simulated_records(i)%trop%ua
        this%simulated_records(i,iens)%trop%va = this%simulated_records(i,iens)%trop%va - this%mean_simulated_records(i)%trop%va
        this%simulated_records(i,iens)%trop%ta = this%simulated_records(i,iens)%trop%ta - this%mean_simulated_records(i)%trop%ta
        this%simulated_records(i,iens)%trop%sh = this%simulated_records(i,iens)%trop%sh - this%mean_simulated_records(i)%trop%sh
      end do
    end do

  end subroutine obs_raob_calc_perb_simulated_records

  subroutine obs_raob_calc_innov(this)

    class(obs_raob_type), intent(inout) :: this

    type(linked_list_iterator_type) it
    integer i, j, k

    ! Innovation = Yo - H mean(X)
    it = linked_list_iterator(this%records)
    j = 0
    do i = 1, this%records%size
      select type (record => it%value)
      type is (raob_record_type)
        call this%innov_records(i)%init_with_size(num_man_level =record%man %num_level, &
                                                  num_sigt_level=record%sigt%num_level, &
                                                  num_sigw_level=record%sigw%num_level, &
                                                  num_trop_level=record%trop%num_level)
        ! Mandatory levels
        this%innov_records(i)%man %p  = record%man %p
        this%innov_records(i)%man %ua = record%man %ua - this%simulated_mean_records(i)%man %ua
        this%innov_records(i)%man %va = record%man %va - this%simulated_mean_records(i)%man %va
        this%innov_records(i)%man %ta = record%man %ta - this%simulated_mean_records(i)%man %ta
        this%innov_records(i)%man %sh = record%man %sh - this%simulated_mean_records(i)%man %sh
        do k = 1, record%man%num_level
          j = j + 1; call check_omb_ua(this%innov_records(i)%man %ua(k), this%mask(j))
          j = j + 1; call check_omb_va(this%innov_records(i)%man %va(k), this%mask(j))
          j = j + 1; call check_omb_ta(this%innov_records(i)%man %ta(k), this%mask(j))
          j = j + 1; call check_omb_sh(this%innov_records(i)%man %sh(k), this%mask(j))
        end do
        ! Significant temperature levels
        this%innov_records(i)%sigt%p  = record%sigt%p
        this%innov_records(i)%sigt%ua = record%sigt%ua - this%simulated_mean_records(i)%sigt%ua
        this%innov_records(i)%sigt%va = record%sigt%va - this%simulated_mean_records(i)%sigt%va
        this%innov_records(i)%sigt%ta = record%sigt%ta - this%simulated_mean_records(i)%sigt%ta
        this%innov_records(i)%sigt%sh = record%sigt%sh - this%simulated_mean_records(i)%sigt%sh
        do k = 1, record%sigt%num_level
          j = j + 1; call check_omb_ua(this%innov_records(i)%sigt%ua(k), this%mask(j))
          j = j + 1; call check_omb_va(this%innov_records(i)%sigt%va(k), this%mask(j))
          j = j + 1; call check_omb_ta(this%innov_records(i)%sigt%ta(k), this%mask(j))
          j = j + 1; call check_omb_sh(this%innov_records(i)%sigt%sh(k), this%mask(j))
        end do
        ! Significant wind levels
        this%innov_records(i)%sigw%p  = record%sigw%p
        this%innov_records(i)%sigw%ua = record%sigw%ua - this%simulated_mean_records(i)%sigw%ua
        this%innov_records(i)%sigw%va = record%sigw%va - this%simulated_mean_records(i)%sigw%va
        this%innov_records(i)%sigw%ta = record%sigw%ta - this%simulated_mean_records(i)%sigw%ta
        this%innov_records(i)%sigw%sh = record%sigw%sh - this%simulated_mean_records(i)%sigw%sh
        do k = 1, record%sigw%num_level
          j = j + 1; call check_omb_ua(this%innov_records(i)%sigw%ua(k), this%mask(j))
          j = j + 1; call check_omb_va(this%innov_records(i)%sigw%va(k), this%mask(j))
          j = j + 1; call check_omb_ta(this%innov_records(i)%sigw%ta(k), this%mask(j))
          j = j + 1; call check_omb_sh(this%innov_records(i)%sigw%sh(k), this%mask(j))
        end do
        ! Tropopause levels
        this%innov_records(i)%trop%p  = record%trop%p
        this%innov_records(i)%trop%ua = record%trop%ua - this%simulated_mean_records(i)%trop%ua
        this%innov_records(i)%trop%va = record%trop%va - this%simulated_mean_records(i)%trop%va
        this%innov_records(i)%trop%ta = record%trop%ta - this%simulated_mean_records(i)%trop%ta
        this%innov_records(i)%trop%sh = record%trop%sh - this%simulated_mean_records(i)%trop%sh
        do k = 1, record%trop%num_level
          j = j + 1; call check_omb_ua(this%innov_records(i)%trop%ua(k), this%mask(j))
          j = j + 1; call check_omb_va(this%innov_records(i)%trop%va(k), this%mask(j))
          j = j + 1; call check_omb_ta(this%innov_records(i)%trop%ta(k), this%mask(j))
          j = j + 1; call check_omb_sh(this%innov_records(i)%trop%sh(k), this%mask(j))
        end do
      end select
      call it%next()
    end do

  end subroutine obs_raob_calc_innov

  subroutine obs_raob_clear(this)

    class(obs_raob_type), intent(inout) :: this

    if (allocated(this%platforms        )) deallocate(this%platforms        )
    if (allocated(this%records          )) deallocate(this%records          )
    if (allocated(this%connectors_man   )) deallocate(this%connectors_man   )
    if (allocated(this%connectors_sigt  )) deallocate(this%connectors_sigt  )
    if (allocated(this%connectors_sigw  )) deallocate(this%connectors_sigw  )
    if (allocated(this%connectors_trop  )) deallocate(this%connectors_trop  )
    if (allocated(this%simulated_mean_records)) deallocate(this%simulated_mean_records)
    if (allocated(this%mean_simulated_records)) deallocate(this%mean_simulated_records)
    if (allocated(this%simulated_records)) deallocate(this%simulated_records)
    if (allocated(this%innov_records    )) deallocate(this%innov_records    )
    if (allocated(this%obs_err          )) deallocate(this%obs_err          )
    if (allocated(this%mask             )) deallocate(this%mask             )

  end subroutine obs_raob_clear

  subroutine obs_raob_init_write(this, key)

    class(obs_raob_type), intent(in) :: this
    character(*), intent(in) :: key

    call fiona_add_dim(key, 'raob_site'      , size=this%platforms%size)
    call fiona_add_dim(key, 'raob_record'    , size=this%records%size)
    call fiona_add_dim(key, 'raob_man_lev'   , size=size(raob_man_p))
    call fiona_add_var(key, 'raob_man_p'     , long_name='raob mandatory pressure'                   , units='hPa'                   , dim_names=['raob_man_lev'                ], data_type='r4')
    call fiona_add_var(key, 'raob_sid'       , long_name='raob station id'                           , units=''                      , dim_names=['sid_len     ', 'raob_site   '], data_type='c')
    call fiona_add_var(key, 'raob_lon'       , long_name='raob station longitude'                    , units='degree_east'           , dim_names=['raob_site   '                ])
    call fiona_add_var(key, 'raob_lat'       , long_name='raob station latitude'                     , units='degree_north'          , dim_names=['raob_site   '                ])
    call fiona_add_var(key, 'raob_elv'       , long_name='raob station elevation'                    , units='m'                     , dim_names=['raob_site   '                ])
    call fiona_add_var(key, 'raob_record_sid', long_name='raob record station id'                    , units=''                      , dim_names=['sid_len     ', 'raob_record '], data_type='c')
    call fiona_add_var(key, 'raob_time'      , long_name='raob observation time'                     , units='hours since 1970-01-01', dim_names=['raob_record '                ], data_type='r8')
    call fiona_add_var(key, 'raob_man_ta'    , long_name='raob temperature on mandatory levels'      , units='degC'                  , dim_names=['raob_man_lev', 'raob_record '], data_type='r4', missing_value=real_missing_value)
    call fiona_add_var(key, 'raob_man_ua'    , long_name='raob wind u component on mandatory levels' , units='m s-1'                 , dim_names=['raob_man_lev', 'raob_record '], data_type='r4', missing_value=real_missing_value)
    call fiona_add_var(key, 'raob_man_va'    , long_name='raob wind v component on mandatory levels' , units='m s-1'                 , dim_names=['raob_man_lev', 'raob_record '], data_type='r4', missing_value=real_missing_value)
    call fiona_add_var(key, 'raob_man_sh'    , long_name='raob specific humidity on mandatory levels', units='m s-1'                 , dim_names=['raob_man_lev', 'raob_record '], data_type='r4', missing_value=real_missing_value)

  end subroutine obs_raob_init_write

  subroutine obs_raob_write(this, key)

    class(obs_raob_type), intent(in) :: this
    character(*), intent(in) :: key

    type(hash_table_iterator_type) site_it
    type(linked_list_iterator_type) record_it
    character(8), allocatable :: cvar(:)
    real(r8), allocatable :: rvar(:), rvar_lev(:,:), rvar_ens(:,:,:)
    integer i, k1, k2, iens

    allocate(cvar(this%platforms%size))

    site_it = hash_table_iterator(this%platforms)
    do i = 1, this%platforms%size
      select type (site => site_it%value)
      type is (raob_station_type)
        cvar(i) = site%name
      end select
      call site_it%next()
    end do
    call fiona_output(key, 'raob_sid', cvar)

    deallocate(cvar)

    allocate(rvar(this%platforms%size))

    site_it = hash_table_iterator(this%platforms)
    do i = 1, this%platforms%size
      select type (site => site_it%value)
      type is (raob_station_type)
        rvar(i) = site%lon
      end select
      call site_it%next()
    end do
    call fiona_output(key, 'raob_lon', rvar)

    site_it = hash_table_iterator(this%platforms)
    do i = 1, this%platforms%size
      select type (site => site_it%value)
      type is (raob_station_type)
        rvar(i) = site%lat
      end select
      call site_it%next()
    end do
    call fiona_output(key, 'raob_lat', rvar)

    site_it = hash_table_iterator(this%platforms)
    do i = 1, this%platforms%size
      select type (site => site_it%value)
      type is (raob_station_type)
        rvar(i) = site%z
      end select
      call site_it%next()
    end do
    call fiona_output(key, 'raob_elv', rvar)

    deallocate(rvar)

    allocate(cvar(this%records%size))

    record_it = linked_list_iterator(this%records)
    do i = 1, this%records%size
      select type (record => record_it%value)
      type is (raob_record_type)
        cvar(i) = record%station%name
      end select
      call record_it%next()
    end do
    call fiona_output(key, 'raob_record_sid', cvar)

    deallocate(cvar)

    allocate(rvar_lev(size(raob_man_p), this%records%size))

    rvar_lev = real_missing_value
    record_it = linked_list_iterator(this%records)
    do i = 1, this%records%size
      select type (record => record_it%value)
      type is (raob_record_type)
        do k1 = 1, size(record%man%p)
          do k2 = 1, size(raob_man_p)
            if (record%man%p(k1) == raob_man_p(k2)) then
              rvar_lev(k2,i) = record%man%ua(k1)
              exit
            end if
          end do
        end do
      end select
      call record_it%next()
    end do
    call fiona_output(key, 'raob_man_ua', rvar_lev)

    rvar_lev = real_missing_value
    record_it = linked_list_iterator(this%records)
    do i = 1, this%records%size
      select type (record => record_it%value)
      type is (raob_record_type)
        do k1 = 1, size(record%man%p)
          do k2 = 1, size(raob_man_p)
            if (record%man%p(k1) == raob_man_p(k2)) then
              rvar_lev(k2,i) = record%man%va(k1)
              exit
            end if
          end do
        end do
      end select
      call record_it%next()
    end do
    call fiona_output(key, 'raob_man_va', rvar_lev)

    rvar_lev = real_missing_value
    record_it = linked_list_iterator(this%records)
    do i = 1, this%records%size
      select type (record => record_it%value)
      type is (raob_record_type)
        do k1 = 1, size(record%man%p)
          do k2 = 1, size(raob_man_p)
            if (record%man%p(k1) == raob_man_p(k2)) then
              rvar_lev(k2,i) = record%man%ta(k1)
              exit
            end if
          end do
        end do
      end select
      call record_it%next()
    end do
    call fiona_output(key, 'raob_man_ta', rvar_lev)

    rvar_lev = real_missing_value
    record_it = linked_list_iterator(this%records)
    do i = 1, this%records%size
      select type (record => record_it%value)
      type is (raob_record_type)
        do k1 = 1, size(record%man%p)
          do k2 = 1, size(raob_man_p)
            if (record%man%p(k1) == raob_man_p(k2)) then
              rvar_lev(k2,i) = record%man%sh(k1)
              exit
            end if
          end do
        end do
      end select
      call record_it%next()
    end do
    call fiona_output(key, 'raob_man_sh', rvar_lev)

    deallocate(rvar_lev)

  end subroutine obs_raob_write

  subroutine obs_raob_final(this)

    type(obs_raob_type), intent(inout) :: this

    call this%clear()

  end subroutine obs_raob_final

end module obs_raob_mod
